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Abstract 


In this report, we have designed an Essentially Non Oscillatory reconstruction for functions 
defined on finite-element type meshes. Two related problems are studied : the interpola- 
tion of possibly unsmooth multivariate functions on arbitrary meshes and the reconstruction 
of a function from its averages in the control volumes surrounding the nodes of the mesh. 
Concerning the first problem, we have studied the behaviour of the highest coefficients of 
two polynomial interpolations of a function that may admit discontinuities of locally reg- 
ular curves : the Lagrange interpolation and an approximation such that the mean of the 
polynomial on any control volume is equal to that of the function to be approximated. This 
enables us to choose the best stencil for the approximation. The choice of the smallest pos- 
sible number of stencils is addressed. Concerning the reconstruction problem, two methods 
have been studied : a first one base on an adaptation of the so-called reconstruction via 
deconvolution method to irregular meshes and a second one that lies on the approximation 
on the mean as defined above. The first method is conservative up to a quadrature formula 
and the second one is exactly conservative. The two methods have the expected order of 
accuracy, but the second one is much less expensive than the first one. 

Some numerical examples are given which demonstrate the efficiency of the reconstruction 
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Dans ce rapport, nous definissons un algorithms de reconstruction essentiellement non 
oscillant pour des fonctions definies sur des mail I ages de type elements finis. Deux probl&mea 
sont Studies : celui de (’interpolation sur un maillage arbitraire d’une fonction (ou I’une de ses 
derivees) pouvant admettre des discontinuity et celui de la reconstruction de cette fonction 
a partir de la don nee de ses valeurs moyennes dans des cellules de controle autour des n 
noeudsdu maillage. 

Deux type d’approximations polynomials sont itudiees : (’interpolation par polynomes 
de Lagrange et une approximation ou Ton impose les valeurs moyennes dans des cellules de 
controle. 

Dans les deux cas, on a etudie le comportement des coefficients du polynome d’appro- 
ximation, en particular dans le cas oit la fonction <711 *.t donnee admet des discontinuity 
situees sur des courbes regulieres par morceaux. 

Ces resultats, qui generalised des ryultats classiques, permettent d’obtenir un critfcre 
de detection des zones de regularity d’une fonction ; dont on a deduit la construction d’une 
interpolation essentiellement non oscillante. 

Le probleme de la reconstruction non oscillante k partir des valeurs moyennes est aussi 
aborde. Deux techniques sont proposees : la premiere qui utilise une methode de deconvolu- 
tion et la seconde qui utilise (’approximation polynomial a partir des valeurs moyennes. Le 
premier type de reconstruction n’est conservatif qu’i, une formule de quadrature pris tandis 
que la seconde Test formellement. Son cout est aussi bien moindre. 

Le probleme du choix d’un ensemble minimal de stencils permettant 1 ’interpolation ou 
la reconstruction est aussi traite. On presente un critere heuristique, qui, pour tous les 
exemples que nous avons traites, semble suffisant. 
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Notations 


• R.[A\ K] : finite dimensional vector space of two variable polynomials over JR, 

• N(n) — ~ ^ : dimension of Ol«(X, K], 

• S ini admissible stencil for solving the Lagrange problem in R„[A', K], see section 1.2, 

• |J£/J| is the euclidian norm of U , 

• D tJ u = where t + j - /, 

• D n u : n*th derivative of u. 

• < u >c 


J u(x)dx 

•* TZkcj 
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Introduction 


During the past few years, a growing interest h\s emerged for building high order accurate 
schemes (i.e of order greater that 2) for compressible flows simulations. It is well known that 
even for smooth initial conditions, these flows may develop discontinuities that make linear 
schemes useless. 

At the beginning of the 80’s, the class of Total Variation Diminishing schemes appeared 
and they have been successfully and widely used with many types of meshes (see for example, 
[1} for a review and, among many others, [2] for simulations on finite element type meshes). 
Nevertheless, one of their main weaknesses is that the order of accuracy falls to first order 
in regions of discontinuity and at extrema, leading to excessive numerical dissipation. 

Various methods have been proposed to overcome this difficulty (adaptation of the mesh 
for example) but one promising way may also be the class of the Essentially Non-Oscillatory 
schemes (E.N.O. for short) introduced by Harten, Osher and others [3, 4, 5, 6, 7). 

The basic idea of E.N.O schemes is to use a Lagrange type interpolation with an adapted 
stencil : when a discontinuity is detected, the procedure looks for the region around this 
discontinuity where the function is the smoothest. Then a reconstruction technique may be 
applied which enables approximation of the function to any desired order of accuracy from 
its averages in control volumes surrounding the mesh points. The approximation is done so 
that it is conservative. 

Some attempts have been made to extend these ideas to multidimensional flows (see for 
example [9]), but only for structured meshes. 

In this report, we intend to study the problem of the reconstruction, up to any order 
of accuracy, of a given function given either by its value at the nodes of a triangulation or 
by its averages on control volumes defined around these nodes so that, in the second case, 
the reconstruction is conservative. This latter problem has already been studied, for smooth 
functions only, by Barth et al. (10]. 

The outline of this report is as follows. In the first part, we consider two ways of approxi- 
mating functions by means of polynomials : the Lagrange interpolation and a approximation 
that conserves the mean on control volumes surrounding given points. In particular, we study 
the problem of the localization of region of smoothness from the coefficients of the polyno- 
mial. We recall the results concerning the approximation of smooth functions by Lagrange 
interpolation and we show that the same results are also true for the second approximation 
we consider in this paper. In the case of unsmooth functions (or more precisely the function 
which regions of possible unsmoothness lies on locally C l curves), we give the asymptotic 
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behaviour of the coefficients of highest degree in the approximation. These latter results 
seem to be new. 

Then we propose an algorithm for LNO interpolation. We try to give some indication 
for selecting the smallest of possible stencils for any order of accuracy. We also propose two 
algorithms for reconstructing a function from its mean values in control volume. The first 
one follows from an adaptation of the so-called reconstruction via deconvolution procedure 
that was originally built for regular meshes. We indicate why, in general, the conservation 
property is formally lost to ensure a high order of approximation when using deconvolution 
technique ; conservation up to a high order quadrature formula is guaranty. Then a secnd 
algorithm, directly inspired from Bath et ai. [10] and Harten et a] (8] is proposed and 
analysed. This latter one ensure formal conservation and, form a computational point of 
view, is much more economical than the first one. 

Some numerical tests are proposed indicate the performance of the various methods. 
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Chapter 1 


Two approximations by polynomials 
in m 2 


1.1 The polynomials in IR 2 


We will denote by IR[X, V'] the vector space of the polynomials of two variables ( X and Y) 
with coefficients in IR. An element of IR[X, K] may be described by its (finite) expansion in 
terms of powers of X and Y : 

/’!■*.>') = £ £ 0. ,X'Y’ (1.1) 

«+jsi «J>0 

The highest integer such that at least one of the coefficients of the monomials X'Y* is non 
zero is called the total degree of P. 

If (x n , »o) >* * point of another expansion of P may be written in terms of the 
monomials (X - io )*(K - yo Y with the help of the Taylor formula. The total degree of P 
does not depend on the point (x 0 , yo). 

In the sequel, we will denote by 1R*(X, V) the (finite) vector space of the polynomials 
of R(X, K) with total degree less or equal to n. This vector space has dimension N{n) = 

( w + l^ w t - D ' a basis of which is the set of monomials (X - Xn) % {Y - yo)* of total degree 
i + ) less or equal to n. 

Let us now describe another interesting basis of R»(X, V). Consider ( A , B,C) a triangle of 
1R J and let us denote by A*, A®, A c the barycentric coordinates of the three points (.A, B, C) 
defined, for any point Af, by : 


M — A a A + As B + Ac C 
A a + A b + A c — 1 


( 1 . 2 ) 


It is easy to see that, for any pair of points, say A and B, the set 
basis of R*[X, >']. 


( A '« A * u » 
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1.2 The approximation problem in 1R 1 2 3 

Throughout this paper, we wilt consider the following approximation problem : let ua con* 
aider the following problem {!) defined for regular enough real valued functiona on II 2 : 

Given N and n, two mtegert, a family of N linear forms defined on 

(^(IR. 2 ) (p > n), find an element P Y] such that for any point i , one 

has L t (P) — "•(**) forl<i<N. 

In the sequel, two kind of linear forma are conaidered : 

• if S {n) = *• 4 family of N points in 1R 2 , the first approximation (Lagrange 

interpolation) is associated with : 

L,(u) - u(A,), 1 <i< N, 


• if M,)i< t <* »» a family of N points in IR 2 and if {<?,, 1 < * < N } ia a fa »ily of control 
volume around the nodes A t , the second approximation we consider is associated with : 

L,(u) =< u > Cl . 

In this second case, we set = (>4 t , C,), <I<JV . This set will be also called a “stencil" 
(though this is not very appropriate in this case) in order to simplify the text. 

In the sequel, we will often make no distinction between i. point element A t of S* n) and its 
coordinates in a suitable frame, (x„ y,). We always assume that the control volume satisfies 
the following properties : 

1. C, f)C, is of empty interior when x ^ j (generally speaking, a collection of segments or 
an empty set), 

2. C, is connex, 

3. we assume that there is a continuous dependency of C , in term of the points of the 

stencil. This guaranties in particular that the area of C'j and the mean values of X*Y* 
on Ci depend continuously of the points of the stencil. 
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2. the following generalized Van der Monde determinant must be non zero : 

Aj... a def \L t (X* Y ')] 

MO M*) MH ••• Li(X') 

L ( 1) L n (X) L„(Y) ••• L(X-) 

We will say that the set is admissible if A$«.» / 0. If a set is admissible, for 

any function u there exist ( n + + V coefficients , (a, ; ), such that the solution of the 

problem (I) is : 

P=i T. +,X'Y>. (..4) 

Is] t ^*1 tj>(l 

In the sequel, the polynomial P will often be denoted as IU(u) in the case of the Lagrange 
interpolation (i.e L,(u) = u(A t )) and fl M (u) in the case of the approximation in the mean 
(i.e L,(u) = < u >c.) 

The problem of characterizing the admissible sets for the Lagrange interpo ation has been 
widely studied, see (11] for example and the references therein. We do not know any %.nteria 
for selecting admissible stencils for the approximation in the mean. 

Remarks : 

1. The condition (1.3) has been jiven for the basis X'Y 1 of F»(X, T]. A similar and 
equivalent condition could have been given for the two other bases we have mentioned. 
The formula (1.4), provided that the monomials X'Y 3 are replaced by the elements of 
the new basis, is also true. 

If card = 3 and if the Lagrange problem is considered, this condition is nothing 
more than the one which says that the three points must not be aligned. 

3. If we were in IR, this determinant, for the Lagrange problem, would be the classical 
van der Monde determinant. 

4. In the case of the Lagrange problem, the set of l n * V -uplets where the con- 

dition (1.3) is not fulfilled is an algebraic curve of IR* 4 " 1 and consequently a closed 
subset of measure zero in R N| *'. In the second case, this i« also true, at least for the 
kind of control volume we consider wich boundary are made of segments, the equation 
of which depend linearly of the coordinates of the nodes of 5 ^ * , . 

In the next section, we address the question of the practical calculation of the coefficient s 


» + ; < n 
1 < l < N 

Li(X n ~ l Y) •••MV'") 
MX^'Y) •••MV'") 
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1.2.1 Determination of the polynomial expansion 

In this section, we use the monomials .Y'K* for expanding polynomials, but any other basis 
would be suitable and the results are immediately transposable. 

The coefficients of the Lagrange expansion are the solution of the linear N(n) x N(n) 
system : 

L t {u) = £ £ a, M x ' yi ) for all 1 <1<N (1.5) 

i=i «+j=i «j>o 

This Lagrange type interpolation problem have been considered by several authors, 
among which Muhlbach [12,13]. 

Since condition (1.3) is true, the Cramer formula applied to (1.5) gives the answer, but 
this solution is not computationally very efficient. Several authors have tried to generalize 
the Newton formula that make the Lagrange interpolation in Bl efficient from a numerical 
point of view, and a very general answer has been given by Muhlbach [12,13]. 

In these papers, he addresses the problem of the “interpolation” by a set of functions 
(/•)•€/ such that the property (X) is true. He calls the set (/*) a Cebysev-system if given 
any function /, for any pair of subsets of /, L and M, having the same (finite) number of 
elements, there exist real numbers a, such that : 

u* = £ a ,£,(/,), for all i € i. (1.6) 

J6M 


For the sake of clarity, we may assume that L — M = {l, • • • N}. He uses the notation : 


h-ft 

Lt 



for denoting the coefficients of f k in the development (1.6). Then, in [13], he shows that if 
one has a Cebysev system (theorem 4.1 pp. 106) : 




fi • • • /•- 1 
Lr-L n 

!]- 

f\ “ * /«-l 
L\- ■ • i 

/] 

Li • ■ L„ 

/ J - 

fl - fn - 1 

Li - ' Ln 

/]- 

A ••/«-. 

L\" • L n - 1 

/] 


(1.7) 


This expression is a direct generalization of the classical Newton formula. Let us make 
sevc. al comments on this formula when applied to our problem : 

1. If one adopts the lexicographic ordering 1, X, Y,X 1 2 ,XY, Y 2 y ■ ■ , X*,X*~ 2 Y, ■ • ,YY W_1 , Y n , 

the previous formula (1.7) must be applied n + 1 times to go from a total degree n to 
a total degree n -f 1. One must also store quite a lot of terms like 


'/»-/* 

L,--L k 
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to build the divided difference table. For example, to go from degree one to degree two, 
one must evaluate and store Cl approximations of gradients by mean of approximations 
on triangles, combine them to obtain approximation on sets of four points (Cl sets), 
five points (C| sets) and six points (one set). Moreover, to go from approximation on 
k points to k + 1 points, (k+ 1) x (k + 1 ^determinants must be evaluated. 


2. From a numerical point of view, the basis X'Y* or (X - a)*(Y — b) 3 are not well suited 
to calculations. This can easily be seen since for any pair (a, 6), 


[£(**’)] 


1 + j < » 
1 < l < N 


[Li\(z - a)* (y - &V]J 


i + > < n 
1 <1 < N 


If (a, 6) is any point of S {r) and if h = mai(„ ,*) € j<«>(|xi - a(, |y ( - 6|), because for 
both kind of linear forms we consider here, we have | Li[(X — a)'(T — 6) J ]| < A ,+ *, 
Hadamard’s inequality shows that : 

|A 5 ,-,| < 

where x(n) = 1 + „ piiLLi^ZLL?) — 0(n*) so that one reaches very quickly 

machine zero though the linear system may be well conditioned. 

An alternative to this last point is to use local coordinates such as the barycentric ones. In the 
E.N.O. method we will develop in a further section, for each point, the natural barycentric 
coordinates are not known a priori so that the work has to be repeated at each interpolation 
call. If this is included in an iterative algorithm, the cost (and the storage) seems to be much 
too important at least for the cases we have considered in this report. For all these reasons, 
we have preferred to use classical inversion techniques for linear systems. 


1.3 A recurrence formula 

In this section, we wish to show another recurrence formula that enables us to obtain the 
coefficients of the expansion of total degree n from those of the expansion of total degree 
less than n in only one step. This recurrence formula may probably be viewed as another 
version of that given in (13], theorem 3.1, page 400 and will be useful in section 1.4. 

Let us begin with some notations. Let Pi, • • •, Ps(n) be a basis of R»(X, Y\ such that 
Pi, • •, Ps( P ), p < n, is a basis of Y). The three basis we have considered in section 
1.2 are of that kind. Let S *** be an admissible stencil. We set : 

Ri = (Li(Pi) ■ • Li(Pn{«))) 

so that the solution of the Lagrange-type problem (I) (where the i^’s are given), 

u,= £ QjLt(Pj) for alii, 
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may be seen as the solution of the linear system 

- Qi/v(«)) = U — (u\ • • ■ «*(„)) , 

where the Z-th row of M is R(. The solutions of system (1.8) are : 

(n) _ det{R\ - -Rj- ■ ■ Rn(*)) 

a ‘ -<ut(Ri - R * •••«#(.>)’ 


where Rt = U denotes the Z-th column. 


( 1 . 8 ) 


(1.9) 


Lemma 1.3.1 Let ( .4, )i <,<*(») be an admissible set in which any of its N(p),p < n, subsets 
is admissible. Then, for a gtven p < n, let / = {i,, • • • **(,,)} and J be ordered sets such that 
1{JJ = {1, • -A r (n)}. If A = (a, ,) is a N(n) x N(n ) matrix, we set 


det(A)i = det(a, ,) 
and det(A)j = det(a t } ) 


1 < i < N(p) ’ 

; e / 

ZV(p) + 1 < i < ZV(n) 
jeJ 


Let us also denote by a\?) the coefficient of P (p) of the Lagrange problem of degree p for nodes 
in I (for degree n, we omit the subscript I). 

Then for any l' < N(p), we have 

Qi** = £/. mr «<(/)=Aro>) Q< f\ 

\ l e - I 1 W f 1 ' 1 ' R* • ’ ‘ Rft(r))i ^(RmpHi ’ ~ ’ Re ■ ■ • Rnm)j ^ 10 ) 

1 1 ' det(Ri • • • ’ 

where o(I) = 1 + p(p + l)/2 + £, € /*\ I n (l 10), Re appears a first time at the l' -th row of 
det(Ri • • • Rf Rn( p ))j and a second time at thel-N(p)-th row of det(R N ^ )+1 •• • R? ■ • RN( n ))j. 


Proof : By switching columns / and l' in (1.9), one gets : 

(«) _ det{R\ • • • Rj ■ ■ Re ■ ■ Rs(n)) 

1 det{R\ • • • Rt • • • Rn(*)) 

Then, a direct application of the generalized Lagrange formula (see [14],pp 19-22) to the 
previous expression gives 

det(/Zj, • • • R l • • • Re • • • Rs («») = 

E/.mr<i(/)=Ar(^j(~ 1 )' T</) det{R\ ■ • ■ Ri ■ • • Rn(j>))i det{Rs(j ,)+ i • ■ • Re • * ■ ffo(»))./ 
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Since any N(p) subset is admissible, one has 

^ _ det(Rj • • • Rt • • • Rn[p))i 
ar detiRx -Rr - Rwn))!' 

and the results follows immediately. • 

Lemma 1.3.2 With the assumptions of lemma 1.8.1, then, for p < n, if N(p) < l < N(n) 
and f < N(p), 

E A ‘/ = o- 

I.c**Kl)=N(p) 

Proof : Apply the Lagrange formula to 

det( R\ • • • Re • • • Rn(?)Rn(p )+ i • * * = 0 » 

and interpret the coefficients in terms of a’s, all equal to 1, and A’s. The lemma 1.3.1 gives 
the result • 


1.4 Approximation of smooth and unsmooth function 
by polynomials 

The problem of interest in this section is the following : Let u be a real function defined on 
an open subset ft of R 2 . We assume that u is n times continuously differentiable on ft except 
perhaps on a subset of ft consisting of a finite collection of locally C 1 curves. Let now T be 
a mesh. For each point of T, we consider one of the two Lagrange-type interpolation of u 
considered above. Is it possible to localize the regions of smoothness of u from the coefficients 
of the polynomial approximation of u ? The answer is yes if additional assumptions are made 
on the mesh. These assumptions guarantee that one can solve the problem ( I) for any order 
from 1 to n, and may be seen as a very natural generalization of classical conditions used in 
the finite elements theory [15]. 

For functions defined on 1R, one knows that the divided differences of u satisfy : 

• If u is smooth on an interval / containing X\, * • • x m , then there exists { € / such that 

/(•>//) 

l*i. *i.- ••.*•!*»] - ■ 


• if u (k) has a jump [u (t) J on /, one has 

[x„x 2 , • • • ,*.|u] = 0((ii w ](*. - *,)*"•). 
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In this section, we intend to generalize these relations, and in particular, the second one 
since this problem seems (surprisingly) not to have been studied yet. The proof appears to 
be rather technical. It is divided into two parts. In the first part, we study the case of a 
stencil >? <n * of <V(n) points where u admits two values, 0 and 1. We give the asymptotic 
behaviour of the highest coefficient of the polynomial of degree n that approximates u when 
it is exactly of total degree n. Then, we define a condition on the stencils that appears to 
be a generalization of the one that says that triangles must not have too small angles to 
ensure a uniform error bound for classical finite elements [15]. Then, using Lemmas 1.3.1 
and 1.3.2, we obtain our result. Let us begin with the case of a stencil in which convex hull 
u is smooth. We first recall the results or Ciarlet and Raviart for Lagrange interpolation and 
we show that their proof is also true for the second kind of approximation, so that a high 
order of approximation is ensured. 

1.4.1 Case of a ’’smooth” stencil 

This problem, for Lagrange and Hermite interpolation, has been studied by for example 
Ciarlet and Raviart in [16]. As we show it in the paragraph 1.4.1, their error study of the 
Lagrange interpolation can easily be extended to the second problem we consider in this 
paper. In the following theorem, K denotes the convex hull of the points of the stencil 
in the case of the lagrange interpolation and the convex hull of the reunion of the control 
volumes C, for the approximation in the mean. 

Theorem 1.4.1 Let <S (n) be an admissible ( for degree n) stencil of IR J , and let h and p be 
respectively the diameter of K and the suprenum of the diameters of the spheres contained 
in K. Let u be a function that admits everywhere in K an + 1th derivative D n * 1 u with 

M n +i = au/>{||Z? n+1 u(x)||;x 6 K) < +oo. 

If n denotes either II or l\u, the projectors on the space of polynomials of degree less than 
n, we have for any integer m with 0 < m < n, 

, U p(||£r U (x) - Z) m n(u)(x)||;x €K)< 

p 


for some constants 

C = C(n,m,S (n) ). 

Moreover, if S^' is obtained from 5 (n) by an affine transformation, then C(n,m,S = 
C(n,m,S< n) '). 

From this inequality, one sees that the "flatter” K is, the poorer the estimation is. A 
direct application of this theorem gives a generalization of our first statement. 
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Proof of the error estimate for approximations in the mean 

The proof of this results is directly inspired form Ciarlet k Flavian without almost any 
change. In all what follows, K denotes the convex hull of Ui<,<v(«) Q. 

Lemma 1.4.2 Let n > 1 be an integer number, let us assume that S (m> = (a,,C<)i <•<*(*) 
is admissible and let u 6 C* with p > n. Then, Jor any x € K and for any integer m, 
0 <m<n, one have : 

D m n{u) = D m u+~ T [D m+1 u(i/,(x))(a. - *)* +1 ] D m Pi (1.11) 

where the p, ’s are the polynomials defined by 

< P. >c,= 

the a, 's are some point of C, and the t/, 's are 


i/,(x) = 8 t x + (1 -0 t )a, 


for 0 < 6 % < 1 

Proof : The proof follows the one given by Ciarlet and Raviart in [16]. 

From the definition of the p,’s, we have : 

nju(u) = < u >c. Pi 

•=i.Ar(«) 

Since C t is connex and u continuous on C„ there exist some a, such that < u >c,— u(a,) 
and consequently : 


N(n) *(»> 

n*(u) = 53 D m n*(u) = 53 u(<*)£ m p. 


•si 


«=1 


Let x £ K. One apply the Taylor formula to order n : 
u(o.) = u(x) + £>u(x).(o, - x) + • • —^( 4 ( 0 , - x) n + ~ ^ D^ufaix)).^ - *)" +1 
where v x satisfies the conditions of lemma 1.4.2. From this, one obtain : 

i>-n M ( U ) = t i £ (o'u(i).(a. - *)'} £ {D“‘ , u(p.(l)).(o 1 - !)•*’} D-p. 

( 1 . 12 ) 
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As in (lo), let us show now that : 


• W<») 

7iE{°'“W(o.-x) , }O m R = 
*• t=l 


0 for 0 < / < m - 1 
D m u(x) for / = m 
0 for m + 1 < / < n 


One first show that 


NM 

L 


is) 



i)'} D~ P , = 0 


for l € {0, 1, • • • ,m - 1}. 


(1.13) 


For m = 0, one notice that 1 = (because FI W let the polynomials of total degree 

less than n > 0 invariant), so that for any constant Ao, 


N{n) 

o = 53 An D m p„ 

•si 

and, in particular, for the constant function A n = u(z). This shows (1.13) for / = 0. 

The proof of (1.13) for / < m - 1 follows by induction : assume (1.13) for / = 0, 1, • • • ,fo 
with lo < m - 2. If A^ is a constant /^multilinear symetric form on R 2 , one can associate 
to it a polynomial P such that D^P = A*, (P = ^AinX* 0 ). From (1.12), and from the 
induction assumption, one has : 


N(n) 

0 *= 53 “ x)* 0 } for allz 6 R 2 

•si 1 

If one sets A ^ one gets (1.13) for 1 = In and the result is true for all 1, 0 < 1 < 

m — 1. 

If / = m, if A m is a constant m-multilinear symetric form on R 2 , since that II u let 
invariant the polynomials of degree less than n, the formula (1.12) combined with what we 
have just shown indicate that : 


N{») 

Am = 53 {^»( a » “ 2 ) < °} 0 m P" for all z € R 2 

• si 

In particular for A m = D m u(x ) and then (1.13) is true for l = m. 

The proof for m H- 1 < n is similar and follows by induction. • 


From that proof, one can see that the key point is the invariance of polynomials by IT W , 
the fact that £* ( " , p l = 1 and the existence of the q,’s. Only that latter point is particular 
to the projection flv we have chosen here. 

The rest of the proof of theorem for n = II M is exactly as in (16], so that we only indicate 
the lemmas that are needed. 
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Lemma 1.4.3 LetSW = {a,, C,} be an admissible stencti LetS {m) = {a*,Cj be an equith 
alent stencil , i.e there exist a matrix B and a point z 0 such that : 

a, s Ba, + *o 

C t = BC, + xo 

Then «S (n) is admissible. 

Lemma 1.4.4 If h is the diameter of K and p is the diameter of the circles contained in 
K , then, if andS^ n ^ are two equivalent stencils, 

lien < | 
iis-'ii < 

P 

The combination of lemmas 1.4.3 and 1.4.4 combined with lemma 1.4.2 gives the theorem 
1.4.1. 


1.4.2 Case of an ’’unsmooth” stencil 

Study of a simplified problem 

Let us consider S (n} an admissible stencil of cardinality L w + ^) &nd So, S\ two non- 

empty subsets of having empty intersection, the union of which is $*"*. Let us consider 
a polynomial P of total degree n such that for all points of So, P has value 0 and for those 
of Si, P has value 1. We define, for e > 0, the set V n of possible stencils for which the total 
degree of P is exactly n : 

Property 1.4.5 The stencil S ( - n> belongs to V' n if and only if, for any vector U which com - 
ponents are either 0 or 1 such that both values are represented, 



R N(n-l) -g- 



(1.14) 


where FU =■ U 

In the definition of property 1.4.5, we have adopted the lexicographic ordering 
Remarks : 

• There are only 2^** -2 such vectors U , so that the number of conditions (1.14) remain 

finite. 
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• The condition (1.14) is homogeneous of degree 0. 

• For first and second order approximation, we have the following result : 

Lemma 1.4.6 If any subset of cardinality (* . + I K A + . V > k < n, of is admissible , 
then the polynomial approximations are of maximal degree for n = 1,2 

Proof : 

- Degree one. The stencil is made of three points A, B,C that make a triangle. P 
is either of type A* or 1 - A* and is of degree exactly one. 

- Degree 2. Let us assume that P is at most of degree 1. Set M = card(Sn) and 
M = S}. We have N + M = N(2) = 6. One may assume that N < M by 
changing P into 1 - P. So, 2N < 6. Then, one can see that there is always, 
for degree 2, at least 3 points that have the same value. Since these three points 
are admissible for degree one and since by assumption, P is either a constant or 
of degree 1 , we see that it must be a constant which is absurd since it takes two 
different values. • 

The same argument applied to degrees n=3,4,5 shows that P must be of degree exactly 
n - 1 because we always have 

(n + 1) (n + 2) n (n + 1) 

4 2 

but fails for degrees greater or equal to 6 (because the previous inequality does not 
hold if n > 5). 

With this in hand, we get the following result, if Ri ; is the following vector : 

R, , = (i,|(X - In )'(K - *«)*] • • £„((* - - tt,)’!) 7 ' . 

Lemma 1.4.7 Let (xn,yn) be an V point of the convex hull K of&*K Assume that the stencil 
satisfies property 1.4.5 for some e > 0. Let P be a polynomial that is 1 on Si and 0 on So. 
If h = max{(|xj - x 0 |, | y t - yo|), (*i.y*) € then there exist two constants Ci(n,e) > 0 

and C 2 (n) > 0 such that the coefficients of the Taylor expansion of P around (x«, y 0 ) • 

P= £ a. ,(*-x„/(K -*«)», 

satisfy 

C 2 (n)h~ n > £ |a. ,!><?,(»»,«)*-. (1.15) 

»+/=n 
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Proof ; We adopt the lexicographic ordering of monomials. The set of admissible stencils 
satisfying property (1.4.5) and |h| < C is a compact subset C of IR**”* because the determi- 
nant is continuous and the values of the mean of X'Y 1 depend continuously of the point of 
the stencil. Let us consider the real functions defined on C by : 

4ij(Au“ *Mn) * J2 \det(Rca • • -U % , • • • 

ial.ll 


and 


V'uMi • • • An) = 


det( Ron • • • V»i 

•Ron) 

det(R y ••• 

Rn) 


in which the vector V %J stands at the "^"-th column and has the value 




and also V t } be £,+,<« P(A %i )U %J where the ”1” is at the L- th position (we refer to the 
lexicographic order). Let us note that P(A t) ) is zero or one, and its value depends only on 
ij and not S (n) . 

It is clear that rp t] is the absolute value of the coefficient of ( X — 2o)*(K — yo Y in the 
Taylor expansion of P. It is also clear that 


/r ( ‘ + '> \det[R« a ---R tr - R 0il ]\<4, r 


by the triangle inequality and because \Li\{X - x n )*(Y - y 0 ) , ]| < h** 1 . So, 


** E M2 

• +JB* 


E 


det(/Znn • • ■ * Kj * • • f^) 

j 


The left hand side of this inequality is a continuous and homogeneous of degree zero function 
on C and hence reaches its minimum. This minimum cannot be zero because that would 
mean that all of the a, ; ’s are zero which is absurd by assumption, and is also independent 
of h (because of the homogeneity property). 

Now let us turn to the second inequality. We have 


E Kil® E IlM* 

ifjs* i+j»* 

The left hand side is also a continuous function on C and is bounded above. Clearly, the 
latter constant does not depend on e. • 


Corollar 1.4.8 With the assumptions of lemma 1.4-1, let n and p be integers satisfying 
p < n. Choose l and V such that N (p - 1) < l' < N(p) and N(n - 1) < l < N(n). Then 
there exist two constants C\ and such that 


C 2 (n)h~ n * 1> > £ \>?, e \>C\(n.€)h-** t 

isardU) =f> 


for any subset of cardinality p. 
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Proof : We apply the definition of A; f and the same techniques as used in the previous 
lemma • 

Now, we may state our main result : 

Theorem 1.4.9 Let S ln) be a stencil satisfying the assumption of Lemma 1.4-7, and for 
some <•> > 0 the property 1.4-5. Let us also assume that the following condition, for some 
0 < a < i 

- >a, (1.16) 

Let u a real function defined on an open subset of ft in H J being C l except perhaps 
on a locally C l curve where its n-th order derivative may have a jump (0"u] such that the 
intersection 1 of that curve and the convex hull of 5 (n) is not empty. Then there exits a 
constant C(n,e 0 , a) > 0 such that the coefficients in the Taylor expansion towards a point of 
2 satisfy : 

£ |a.,|>C(n 1 c n ,a)^. (1.17) 


Remark : In the case of the Lagrange interpolation, for n=l, the condition (1.16) just says 
that the smallest angle of the triangle is not too small. 

Proof : 

Let u be a real valued function that admits continuous derivative up to order p — 1 and a 
jump in its p th derivative across a locally C l curve. C. Let us remind the recurrence relation 
of lemma 1.3.1 : 


,<"> - 


= £ 


l.rnr4(l) mp 


K l tt i7 


where the subscript / (reap. () corresponds to (i, j) with i + j = n (resp. * + j = p) in the 
lexicographic ordering. 

For the sake of simplicity, we assume that the curve C divides K in two parts, K * and 
K~ , each of them being on one side of C. Let (lo.Vo) he a point of K, say in K*. We 
establish first the proof of the theorem for the Lagrange interpolation and show at the end 
of it how to adapt it for the approximation in the me^n. The Taylor formula gives, for 


“(a) = EfJ ^ D % ,u{x a ,ya)(ar - x 0 )'(s “ VoY 


+ pT A,u(t/(X n, y - 0))(a, - Xo)‘(a* - y o y 

where ^(x 0 , yn) stands for (i/ xn+( 1 -u)a t , u yo+(l for some v € JO, 1], D t) u(u(x Q ,y 0 )) — » 
For a point of $ ( “ , n # + . one gets : 

r-i j j 

u (°) = Y. 7i £ D tJ u(x n ,ya)(a, - x^o* - y„Y + - ]T A '>( a ’ ~ ~ yof (118) 

lud ' |+J»1 P‘ !♦;«(! 
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where A %} —* D~. In order to get this result, on has to consider, on the straight line 
[(x 0 , y 0 ), al, the function v defined by : 

• if x € |(x„.y 0 ),a] is between (in,Jfo) and x c = [(*n,yn).a] f)C\ v * u. 


• if x € |(xo,yo),a] is between x c and a, 

w = « + A £ ( A,](x - x„)’(y - y n y 
P- 


where [£>,,] represents the jump of D t) u at xc- 

This function admits a continuous p tk derivative over [(in, J/n).a], so one obtains the result 
1.18. 

Let ( > 0 be given. If the diameter of K is small enough, from these relations and the 
rules of evaluation of determinan ts, one gets : 



d:,u\(x -*»)•(/ -!»)'] 

D*L_[(X 

/Ion D t} Lm+ 1 [(•' “ *o)'(y “ Ito) 1 ) “■ ^<>p 

D~L N(r) [(X - xqY(Y - yoY) 
det(Roo • • - Rap) 


(1.19) 



where the index l' stands for («,;) in the lexicographic order, L,(u) = u(a.) and |o(l)| < «. 

For the sake of simplicity, let us denote by m and >// the coefficients obtained from the 
first part of the right hand side of equation 1.19 by replacing D * by 1 and 0" by 0 for m 
and vice versa for ui so that 

= (oj +«(!))/*/ + (b; + «(i))»i. 

We have to notice that the sum of m and v/ is one, and that if all the points of / are on the 
same side of C, then either m or is zero. Then, using lemma 1.3.2 and the above remarks, 
we have : 

«)*' = [£>.,] Z + Z <*> )/•/*'/+ Z oOwa 1 /, 

rar4(l)xf i rariUDxp rar4(l)mp 
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to that 


114,11 V~ 


E 


i//A/ 


< A r “*|a<| + € k*~* £ 

ntr4(f)t 


(Im/I + M)I A / 


i ft 


(1.20) 


// 


Because of inequality (1.16) and front Hadamard's inequality, one has : 


Im/| + K< < 



- **Y(Y - + / £ £.((* - *«)'(y - *)f 

V«€X- 


(Note that we have identify the point of the point a and its index in the ordering.) This 
latter expression is bounded above by a constant C n because F corresponds to (i,j) in the 
lexicographic order, without any additional conditions on the geometry and because L t is 
continuous. With the help of lemma 1.3.2, the expression // of (1.20) is bounded above by 
a constant C\ . 

The term / of (1.20), with the help of lemma 1.3.1 can be seen as the P h coefficient of 
a polynomial that admits the values L m \(X - in)'^ - Vo)*] where a 6 K~. At least one of 
these value is not zero because the curve C cuts K so that there are points of on both 
sides of C and (x 0 ,yo) € K*. When summing up these equalities for i + j = n, the same 
proof as in lemma 1.3.1 shows that this expression is bounded above by a constant C; / 0. 
Then, one gets our result for < small enough. 

In the case of the approximation of the mean, the same kind of arguments may be applied 
because : 


• if C A C K* or A", these exist a,6C t such that 


< u >c,= t*(a.) 


• if C„ cuts C, then there exists a+ € C 4 f) K* and a" € C m f| K~ such that : 


< u > c ,= 


aire{C u (\K*) 

aire(C 4 ) 


«*(<*•) + 


a\re 


(C.f\K 


aire(C 4 ) 


-«(*:)• 


This result enables us to detect the regions of smoothness from those where a jump in 
one of the derivatives occurs. 
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Chapter 2 

An E.N.O. Reconstruction Technique 


In recent papers. Harten and several other authors (3,4 ,5,6,7] have tried to derive numerical 
methods that are able to achieve a higher order of accuracy than classical TVD methods. 
The.e are several versions of these techniques, but they can be generally viewed in the 
following way : starting from some approximation of a real function u (point values or 
average values in some control volumes), find a pointwise high order approximation u. Two 
tools are then used : 

• an essentially non oscillatory Lagrange interpolation of a function tv, 

• this function tv may be u itself if one starts from point values or either the primitive 
function of u or its convolution product with the characteristic function of a copy of 
the control volume if one can pass from one to another by a translation 

The latter deconvolution technique can only be applied, at least in its standard version, to 
regular meshes as shown in section 2 . 2 . 1 . 

In this section, we want to adapt both tools to the situation of unstructured meshes. At 
least for the second point, the situation seems at first glance very bad : the reconstruction 
via primitive function cannot be applied in the case of unstructured meshes because the only 
solution would be to apply it to integrals over domains like Vy = [aj, 6 j] x [ 02 , 63 ], possibly 
with a suitable transformation of the plane as in [9], which are not in general the union of 
control volumes. 

Now, the reconstruction via deconvolution technique can only be applied to regular 
meshes (i.e. meshes where the control volumes are translated form one node to another). In 
this section, we show how to adapt this technique for irregular meshes. 

In what follows, T is a triangulation of Q, a domain of R 2 , u is a function defined on 
that domain. Around each node t of T, we may define a control volume C, in many different 
ways. An example is (see Figure 2.2) the control volume whose boundary is the segment 
joining the centroids G of the triangles ( i,j,k ) having i as a vertex and the middles /j , /* 
of the segments of those triangles (type I). Another type of control volume is obtained by 
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considering each triangle as the control volume of its centroid. The triangulation we need 
to describe the ENO algorithm is not T but another one built from the centroids of the 
triangles of T. 


2.1 Non oscillatory interpolation 

Let us consider n > 0. In this section, we show how to generalize the n + 1-th order E.N.O. 
interpolation technique exposed in [4,5] to unstructured meshes. As the results of chapter 
1 indicate it, we must deal with meshes where the stencils we need satisfy the property of 
theorems 1.4.1 and 1.4.9. This will be the case for most meshes. 

Let T be a triangulation of ft, a domain in H 2 , and u a function defined on that domain- 
The results we have obtained in chapter 1 can be summarized as follows : 

• if «S (B * is an admissible stencil such that u is smooth in its convex hull, then 

E K ,i. 

•+j=# 


remains finite, 

• if in the convex hull of S (n \ u admits continuously differentiable derivatives only up 
to the order k < n, then 

£ k,l = o([u<->] A*-), 

«+/=« 

where h is the diameter of S^ n K 

Then, as suggested by (4,5], the E.N.O. algorithm we propose is to consider IIi(u) defined 
on C, by the following recursive algorithm : 

For a node i, 

1. Let { Tj } be the set of triangles of T having * as a vertex. Consider all the linear 
interpolations where the T t 's are the stencils. Choose the one, T mtn , where the sum 

E K it. 

«+ ; =i 

is minimal. We set «S (1) = the nodes of T^. 

2. Let S* n ~ l * be the stencils defined at the previous step. Consider all the nodes sur- 
rounding «S (o ~ 1, in T and consider all the stencils obtained form by adding n+ 1 

of the nodes surrounding S*"' 1 '. Choose the stencil minimizing : 

L M* 
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We have intentionally left the second point imprecise because it is obvious that the 
number of stencils to consider is in general huge. To give an example, if one considers Figure 
2.3 for which n = 2, one sees that possible stencils are the vertices of triangles 7 L and 3 of 
the ten other triangles. This can be repeated for each of the three edges of 7 'm»n and leads 
to a total of 3 x C\ a — 360 possible stencils ! So, one has to define criteria for choosing the 
"good” and "bad” stencils. These criteria are essentially heuristic and a priori ones . 

One that seems natural is that when one considers the control volume around each node, 
the collection of the control volumes of all points of the stencil should be convex. Another 
one is that the criteria leads to the smallest possible number of stencils, but the stencils must 
not be confined in a particular angular area of the plane, in order not to favor any direction. 

With this in mind, two possible sets of stencils for third order interpolation are (see 
Figure 2.3) : 

• the nodes of triangle plus, for each of its edges, the three additional nodes of 
triangles Ti, 72, T 3 . This leads to a maximum of three stencils per triangle, 

• or the nodes of triangles of T mtn plus, for each of its edges, the three additional nodes 
of triangles of 

- t„t 2 ,t 3 , 

- T \ , T 3 and T 4 or Tj, 

- Ti,T 3 and T 6 or 7V, 

- Ti, T 0 or Tio and one of the six triangles Ti,T 3 ,T T 5 , T 6 ,Tt- 

The second solution leads to a maximum number of 52 stencils once T«u, has been found. 
We have made several tests to evaluate the "performance” of each type of stencil. They are 
given in section 3 and indicate that the first kind of stencil is sufficient. 

These experimental results suggests the following algorithm for selecting stencils for 
higher order of accuracy (of course, in this case, the problem of choosing the smallest number 
of possible stencils is much more critical than for third order interpolation !). 

To begin with, we notice that the six points of T mn , Tj, T 2 and Tj, denoted by in the 
sequel, form a set isomorphic to the standard stencil for P 2 interpolation (see Figure 2.4). 
Now, this stencil can be embedded in the standard stencil for P 3 interpolation (Figure 2.4) 
where the additional points 7,8,9,10 are the new points introduced for 4‘* approximation. 

Let us consider a “side” of say the points 4,5,6. Since the mesh is conformal, the 
segment 4-5 is the edge of two triangles in general, T\ and a new one, T 4 . The same is 
true for the segment 5-6 and a new triangle T% is introduced. Then the sides of these new 
triangles are the edge of some triangle, T6 and possibly T% for T 4 . So, at most four new 
triangles, T 6 , T 7 , T$. T 0 are introduced and the new stencils are made of and the points 
of {T 4 ,T 6 ,T 8 ,Ts} or T 4 ,T 5 ,T 0 ,T 7 } or {T 4 ,T 8 ,T,,T 5 } when T 8 # T 0 These two or three 
new stencils are isomorphic to the standard P 3 stencil as shown on Figure 2.4. These new 
procedure may be repeated for the three sides of and so at most 12 stencils must be 
considered for «S (3) . 
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" y ?- 


By using the similarity between «S < " ) and the standard stencil for P n interpolation, in- 
duction can easily be applied to construct a small, but probably sufficient, set of stencils 
S< n + V 

This particular interpolation is n -f 1-th order accurate; because u is a polynomial of 
degree less than n, we have fli(u) = u. This property ensures the n + 1-th order accuracy 
{16}, and in particular, we have the estimations of theorem 1.4.1. 


2.2 Two conservative ENO reconstruction 

Two kinds of reconstructions are now considered. A first one inspired from the so-called 
reconstruction technique and a second one, already considered by Harten et al. without 
rigorous justification in a very recent paper [8]. 


2.2.1 Deconvolution technique revisited 

If (x), e K is a regular mesh of R and u is a real valued function on IR, the reconstruction by 
the deconvolution technique consists of applying the previous algorithm, not to u but to 

ti.(y) = / u(x + y-x t )dx, (2.1) 

*•+1 x i ■'*»-!/» 


where, as usual, the mesh size Ax = Xi +i - Xi is constant and x, + j/ 2 = x< + Ax/2. In 
particular, we see that u, does not depend on t and that u(x,) is the average of u on 
(x,-i/ 2 ,x, + i/ 2 ]- These values are assumed to be known. Let FI j (tZ) be the m + 1-th order 
Lagrange interpolation as described in the previous section, with m > n. Then, the idea is 
to perform a Taylor expansion of u and its successive derivatives around x,, to truncate them 
at order n - k, to replace the values E, • • •, E* B> by IIi(E) and its n successive derivatives, 
and to replace the values of u by those of n 2 (u), the approximation of ti we are looking for : 


n.(s)(ft) = Sr a ,.»‘ n - ^rH 





( 2 - 2 ) 


where 


/’■*'"( X - x.)‘ dz 

J 


a. = 


Ax 


The linear system is easily invertible because the matrix is upper triangular and its diagonal 
consists of all “l”s. Furthermore, it is shown in [4], for example, that the average value of 
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n 2 (u) over [*i-j/ 2 , 1,-h/j] is exactly E,. Last, this approximation has the desired order of 
accuracy when u is smooth because polynomials are left invariant by the construction. 

This latter point is the fundamental reason why one achieves the expected order of 
accuracy . Polynomials are left invariant by the construction because the shape of the 
control cells does not change from one point to another. If this where not the case, i.e if 
Ax , +j /2 = x.+j - x, were not constant, the formula (2.1) would indeed depend on the point 
x, and this property would be lost. In order to show this, we simply consider u(x) = x and 
a m + 1-th order interpolation that has values u* at points x t . Assuming that {xo, Xj, • * *} is 
the stencil selected by the E.N.O. algorithm, we have : 

Ax 3 /2 + Ax.1/2 


IIj(E) = En + K(x — Xjj) + • • • where 




Ax 


1/2 


When the mesh is not regular, K yt 1 in general. To obtain Il 2 (u) = u, one must have : 

n,(u)(y) = u(y) + a,u'(y), 

n 2 (u/(y) = K = Ay)- 

The second equation indicates that one must have K = 1 which is, in general, not true. 

To overcome this problem, we propose the following technique : apply the ENO search 
algorithm not to u defined by equation 2.1 but to : 


E(y) = 


area(C s t 


Cymt) b s(n) 


,(x + y - x 0 )dx, 


(2.3) 


where is any possible stencil around the node Xo that one has to test and C si ») is the 
union of the control volumes of each node in «S* n> : 


C S ymj — (J Cj. 


(2.4) 


Now, one has to evaluate the integral (2.3) from the average value of u. This cannot generally 
be achieved for any function, but is possible for the polynomials of IR^X, V], at least in 
general. This will be true for the N(n) linear forms over R»(A\ Y\ when 


< P > c = — 1 — f P( x ) dx % for all P 6 K], 

areau, Jc, 


(2.5) 


are independent. For all the meshes we have considered, these linear forms where always 
independent, so the problem had a solution. If this is true, then one can find coefficients 
a i(y)> 1 < f < ^( n ) 80 that 


1 / N{n) 

— — r / u(x + y - xo) dy = £ a t (y) < u > Cl , 

area(C s i«») Jc s(m> 


( 2 . 6 ) 


when u belongs to R,[X,y). If not, the equation (2.6) gives a n + 1-th order quadrature 
formula. With all this, we get the following theorem whose proof is obvious. 
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Theorem 2.2.1 The algorithm defined by the E.N.O. technique with (2.3), (2.4), (2. 5) and 
(2.6) leaves invariant the polynomials of 1R 1I [A',K] and hence gives a n + 1-th order ap- 
proximation of smooth functions. Moreover, this approximation is conservative up to the 
quadrature formula (2.6). 

2.2.2 A formally conservative reconstruction 

A second reconstruction technique, which is formally conservative, may be the following : as 
noticed in theorems 1.4.1 and 1.4.9, the regions of smoothness of afunction u may be detected 
from the behaviour of the highest order coefficients of the approximation in the mean I1 a/ ( u) 
for all the possible stencils surrounding a given node s. The reconstruction we consider is 
the defined by the polynomial having the lower highest order coefficients because, as in the 
previous section, the results of chapter may be summarised as follow for the approximation 
in the mean : 

• if S is an admissible stencil such that u is smooth in its convex hull, then 

51 k/|. 

i+j=n 


remains finite, 

• if in the convex hull of S (n \ u admits continuously differentiable derivatives only up 
to the order k < n, then 

«+>=« 

where h is the diameter of S (n) . 

Then exactly the same algorithm as before may be applied in order to select the “smoothest” 
stencil for the approximation around a given node. Since this node always belongs to any 
possil ie stencil, then the reconstruction is conservative from the definition of I1 a/. For the 
same reasons than in the previous paragraph, the same accuracy properties are also true, 
and the constant of theorem 1.4.1 is the suprenum of all the 6'(n,«S (n) ) for all the stencils 
considered in the reconstruction. 

In practice, the set of possible stencils is defined by induction from the algorithm we 
propose in paragraph 2.1. 

2.3 Some remarks for the practical calculation of the 
reconstruction 

In section 1.2.1, we have discussed the problem of the practical determination of the coeffi- 
cients of a Lagrange interpolant because the linear systems to be solved have in general small 
coefficients. The same problem also arises here if one uses the monomials ( X -xo )'(Y -yoY to 
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determine the ft|(y) of equation 2.6. This can easily be seen by using Hadamard's inequality 
as in section 1.2.1. 

For a given node x n , the E.N.O. technique we propose naturally introduces one triangle 
having that node as vertex, the triangle T mn as in Figure 2.3 . So, as in section 1.2.1, we 
will use the barycentric coordinates towards that triangle for practical computations. 
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Chapter 3 

Numerical examples 


We have performed several tests on the second, third and fourth order E.N.O. interpolation 
and E.N.O. reconstruction, but we only report the third and fourth order results since they 
are a priori more challenging. In particular, we intend to check numerically that the expected 
order of accuracy is in fact reached for smooth functions. 

In all these examples, we have assumed that the control volumes are of "type I”. The 
practical calculations of the averages in these control volumes have been performed with a 
5-th order quadrature formula (15]. 

The tests on smooth functions are performed on : 

u(x, y) = cos( 2ir(x 2 + y 2 )). 


The mesh size h has been measured by choosing the largest segment of the tri angulation. 
All the error estimates have been obtained on irregular meshes as the one presented on Figure 
3.5. These meshes are obtained as random perturbations of regular structured meshes. The 
set of points that one obtains is triangulated by the Brower algorithm to get a Delaunay 
triangulation. The main difference between such a mesh and the regular structured one is 
that the number of triangles each node belongs to is different. We also have done the same 
tests with regular meshes, and we have not seen any degradation of the convergence. 

The locally smooth function we have chosen is obtained by a modification of that used 
by Harten in [7] for example : if <p is any angle, let be : 


/♦(*.]/) 


if r < f 4 (x,y) - -rsin(|r J ), 

1 if r > i, f 4 (x y y) = 2r - sin (3*r), where r = x + tan (^)y, 

. ‘f l r l<5> /♦(*» y) = I sin (2wr)|, 


and let u be : 


(3.1) 


{ if x < j cos xy, 
if x > ^ cosxy, 


**(*•*) aB / v G75(*‘¥)' 

*(*. y ) = y) + c °s (2 xy). 


(3.2) 
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The function defined by (3.1)-(3.2) shows discontinuities in the function itself and its first 
order derivatives; some of the discontinuities are straight lines (never aligned to the mesh), 
one is a curved line where the jump changes from one point to another. Last, the behaviour of 
u is basically one-dimensional on the left of the curve x = cos icy/2 and really two-dimensional 
on the right. 

A plot of this function is given in Figure 3.6. One should obtain straight lines and smooth 
discontinuity transitions contrary to what is shown in the Figure : this is an effect of the 
graphic device adapted to Pi interpolation. 

Third order interpolation and reconstruction No visible difference can be seen be- 
tween the two kind of reconstructions (directly from mean values and by deconvolution), 
except on the cost of each algorithm. Because of that, most of the results are presented from 
the “reconstruction on the mean” technique except if specified. 

The two types of stencils (as presented in section 2.1 ) have been tested in this case. The 
use of the second type of possible stencils results in a much more expensive approximation 
(in general, one must test 52 stencils per triangle versus only 3 in the simplest version) and 
the results have never been dramatically improved. All the results that are presented bellow 
have been obtained with the 3 stencil version of the method. 

We have displayed in Figure 3.7 the L* error of the interpolation for the smooth test 
case. The plain curve with squares is obtained with the E.N.O. interpolation, the plain curve 
with circles is obtained with the E.N.O. reconstruction. The dashed line indicates the slope 
—3. One can see that the expected order of accuracy is indeed reached. 

In Figures 3.9-(a) and 3. 10-(a), we have displayed the node values of the E.N.O. recon- 
struction for two meshes (1600 nodes and 6400 nodes). To better see the behaviour of both 
approximation techniques, we also present cross-section on three lines : Y = 0.75, Y = 0 and 
Y = -0.45 (Figures 3.12, 3.13 and 3.14. The approximations are obtained from the 1600 
nodes mesh (see Figure 3.5) by the deconvolution technique. The latter line goes through 
one of the triple points (see Figure 3.6). One can see that the various discontinuities and 
the smooth regions are well captured by both techniques. 

The Figure 3.11 represents the overall values of the reconstructed function by deconvo- 
lution and can be compared to Figure 3.10-(a). 

Fourth order interpolation and reconstruction The same tests have been performed 
for the fourth order interpolation and reconstruction of the smooth function. The Figure 3.8 
shows the L x error of the ENO reconstruction for random meshes. The Figures 3.9-(a) and 
3.10-(b) shows an overall picture of the fourth order ENO reconstruction. A better capture 
of the area surrounding the triple points is the only visible difference between third and 
fourth order reconstruction. 


To end this section, we must note that the algorithm for choosing the stencils may lead 
to some difficulties at the boundaries as can be seen in Figure 3.6 on the left upper corner : 
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the most left upper triangle of the mesh (Figure 3.5) does not admit any additional points 
of the type we consider to make a stencil. 


Conclusion 


In this report, we have developed and analyzed three methods for the reconstruction of a 
function admitting discontinuities only on regular planar curves from their node values or 
from their averages in control volumes that surround them. They are based on two polyno- 
mial approximations of real valued functions, the Lagrange interpolation and what we have 
called the approximation in the mean. In order to give a firm basis to the Essentially Non- 
Oscillatory interpolation and reconstruction, we have studied the behaviour of the highest 
order coefficients of these two polynomial approximations of smooth functions and unsmooth 
ones for which the discontinuities lie on regular curves. We have also given an adaptation 
of the so called ’’reconstruction via deconvolution” to irregular triangulated meshes. The 
approximation in the mean leads to a very efficient ENO reconstruction that seems to be 
the cheapest one. We also propose an heuristic algorithm for selecting the smallest number 
of possible stencils. 

These techniques have been shown to work quite well on smooth and unsmooth functions. 
In particular, we have shown in these examples that the minimum number of possible stencils 
was sufficient for our purpose. We have not been able to make a rigorous study of the quality 
in term of convergence of these new approximation; that study is clearly needed. In a future 
paper, we will give application of these new approximation for deriving new algorithm for 
compressible Fluid Mechanics. 
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Figure 2.3: Some possible interpolation points. Circles : points of T** (second order in- 
terpolation), black circles : points that may be added to obtain a stencil for third order 
interpolation. 



Figure 2.4: Construction of a stencil for fourth order approximation : the black circles 
represent the points to be added to a stencil for third order approximation (circles). The 
standard P 3 stencil is shown. 
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Figure 3.5: Typical mesh. 1600 nodes, 3042 triangles. 
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Figure 3.6: Exact function. Min=-1.331, Max=2.650, 6 = 0.1. 
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Figure 3.7: L 00 error for f(x,y) = cos (27r(i 2 + y 2 )]. Squares : E.N.O. interpolation only, 
Circles : E.N.O. + reconstruction. Dashed line : slope -3 
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0 



Figure 3.8: L® error for /(x,y) = cos [2x(x J + y 2 )]. Squares : Dashed line : slope -4, 
D ashed -f squares : E.N.O. + reconstruction 
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(a) T hird order 




(b) Fourth order 


Figure 3.9: E.N.O reconstruction “in the mean” with the 1600 nodes mesh. Min = -1 30 
Max=2.6, 6 = 0.1. 43 



(a) Third order 



Figure 3.10: E.N.O. reconstruction “in the mean’ 1 with the 6400 nodes mesh. Min=-1.30, 
Max— 2 6, 6 = 0.1. 44 





Figure 3.11: 3th order E.N.O. reconstruction with the 6400 nodes mesh from the deconvo- 
lution technique. Min=-1.325, Max=2.650. 6 = 0.8281 10 -J . 
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